COMMENT ⊗   VALID 00010 PAGES
C REC  PAGE   DESCRIPTION
C00001 00001
C00002 00002	% This program by D. E. Knuth is not copyrighted and can be used freely.
C00004 00003	@* Introduction.
C00011 00004	@* The character set.
C00016 00005	@* Inputting the data.
C00019 00006	@* Approximating the darkness levels.
C00026 00007	@* Error diffusion.
C00031 00008	@* Computing the diffusion tables.
C00037 00009	@* The main program.
C00038 00010	@* Index.
C00049 ENDMK
C⊗;
% This program by D. E. Knuth is not copyrighted and can be used freely.

% Here is TeX material that gets inserted after \input webmac
\def\title{BTONE}
\font\logo=logo10
\def\MF{\logo METAFONT}
\magnify{\magstep1}
%\pagewidth=4.2truein % estimate to match CACM line length
\setpage
%\tolerance=1000

%\advance\topskip by \baselineskip	% doublespacing
%\advance\smallskipamount by \baselineskip
%\advance\baselineskip by \baselineskip

\def\con{\par\vfill\eject % finish the section names
  \rightskip 0pt \hyphenpenalty 50 \tolerance 200
  \setpage
  \output{\normaloutput\page\lheader\rheader}
  \titletrue % prepare to output the table of contents
  \pageno=\contentspagenumber \def\rhead{TABLE OF CONTENTS}
  \message{Table of contents:}
  \topofcontents
  \line{{\bf Sample}\hfil Section}
  \def\Z##1##2##3{\line{\ignorespaces##1
    \leaders\hbox to .5em{.\hfil}\hfil\hbox to2em{\hss##2}}}
  \readcontents\relax % read the contents info
  \botofcontents \end} % print the contents page(s) and terminate
@* Introduction.
This program prepares a \MF\ program for a special-purpose font
that will approximate a given picture, one ``character'' per line.
The input is assumed to have the form described in the \\{halftone}
program. The output will say `\.{input btone}' followed by a
sequence of lines like
$$\hbox{\tt row(10); cols(3,15,16,17);}$$
this means bits 3, 15, 16, and 17 of the character for row 10 should be black.

@ Here's an outline of the entire Pascal program:

@p program btone(@!input,@!output);
label @<Labels in the outer block@>@/
const @<Constants in the outer block@>@/
type @<Types in the outer block@>@/
var@?@<Global variables@>@/
@#
procedure initialize; {this procedure gets things started properly}
	var@?@<Local variables for initialization@>@/
	begin @<Set initial values@>@;
	end;@#
@<Basic procedures@>
begin initialize; @<The main program@>;
end.

@ The picture in the input data must contain fewer than |max_m| rows and
|max_n| columns.

@<Constants in the outer block@>=
@!max_m=200; {$m$ should be less than this}
@!max_n=200; {$n$ should be less than this}
@!big_n=1600; {$8n$ should be less than this}

@ The main program has one statement label, namely |cleanup_and_terminate|.

@d cleanup_and_terminate=9998
@d finish==goto cleanup_and_terminate
 {do this when all the pictures have been output}

@<Labels in...@>=cleanup_and_terminate;

@ It's convenient to declare a macro for incrementation.

@d incr(#) == #←#+1
@* The character set.
We need translation tables between ASCII and the actual character
set, in order to make this program portable. The standard conventions of
{\sl \TeX: The Program\/} are copied here, essentially verbatim.

@d text_char == char {the data type of characters in text files}
@d first_text_char=0 {ordinal number of the smallest element of |text_char|}
@d last_text_char=127 {ordinal number of the largest element of |text_char|}

@<Types...@>=
@!ASCII_code=0..127; {seven-bit numbers}

@ @<Glob...@>=
@!xord: array [text_char] of ASCII_code;
	{specifies conversion of input characters}
@!xchr: array [ASCII_code] of text_char;
	{specifies conversion of output characters}

@ @<Set init...@>=
xchr[@'40]←' ';
xchr[@'41]←'!';
xchr[@'42]←'"';
xchr[@'43]←'#';
xchr[@'44]←'$';
xchr[@'45]←'%';
xchr[@'46]←'&';
xchr[@'47]←'''';@/
xchr[@'50]←'(';
xchr[@'51]←')';
xchr[@'52]←'*';
xchr[@'53]←'+';
xchr[@'54]←',';
xchr[@'55]←'-';
xchr[@'56]←'.';
xchr[@'57]←'/';@/
xchr[@'60]←'0';
xchr[@'61]←'1';
xchr[@'62]←'2';
xchr[@'63]←'3';
xchr[@'64]←'4';
xchr[@'65]←'5';
xchr[@'66]←'6';
xchr[@'67]←'7';@/
xchr[@'70]←'8';
xchr[@'71]←'9';
xchr[@'72]←':';
xchr[@'73]←';';
xchr[@'74]←'<';
xchr[@'75]←'=';
xchr[@'76]←'>';
xchr[@'77]←'?';@/
xchr[@'100]←'@@';
xchr[@'101]←'A';
xchr[@'102]←'B';
xchr[@'103]←'C';
xchr[@'104]←'D';
xchr[@'105]←'E';
xchr[@'106]←'F';
xchr[@'107]←'G';@/
xchr[@'110]←'H';
xchr[@'111]←'I';
xchr[@'112]←'J';
xchr[@'113]←'K';
xchr[@'114]←'L';
xchr[@'115]←'M';
xchr[@'116]←'N';
xchr[@'117]←'O';@/
xchr[@'120]←'P';
xchr[@'121]←'Q';
xchr[@'122]←'R';
xchr[@'123]←'S';
xchr[@'124]←'T';
xchr[@'125]←'U';
xchr[@'126]←'V';
xchr[@'127]←'W';@/
xchr[@'130]←'X';
xchr[@'131]←'Y';
xchr[@'132]←'Z';
xchr[@'133]←'[';
xchr[@'134]←'\';
xchr[@'135]←']';
xchr[@'136]←'↑';
xchr[@'137]←'_';@/
xchr[@'140]←'`';
xchr[@'141]←'a';
xchr[@'142]←'b';
xchr[@'143]←'c';
xchr[@'144]←'d';
xchr[@'145]←'e';
xchr[@'146]←'f';
xchr[@'147]←'g';@/
xchr[@'150]←'h';
xchr[@'151]←'i';
xchr[@'152]←'j';
xchr[@'153]←'k';
xchr[@'154]←'l';
xchr[@'155]←'m';
xchr[@'156]←'n';
xchr[@'157]←'o';@/
xchr[@'160]←'p';
xchr[@'161]←'q';
xchr[@'162]←'r';
xchr[@'163]←'s';
xchr[@'164]←'t';
xchr[@'165]←'u';
xchr[@'166]←'v';
xchr[@'167]←'w';@/
xchr[@'170]←'x';
xchr[@'171]←'y';
xchr[@'172]←'z';
xchr[@'173]←'{';
xchr[@'174]←'|';
xchr[@'175]←'}';
xchr[@'176]←'~';@/
xchr[0]←' '; xchr[@'177]←' ';
	{ASCII codes 0 and |@'177| do not appear in text}

@ @<Local variables for init...@>=
i:0..last_text_char;

@ @<Set init...@>=
for i←1 to @'37 do xchr[i]←' ';
for i←first_text_char to last_text_char do xord[chr(i)]←@'177;
for i←1 to @'176 do xord[xchr[i]]←i;
@* Inputting the data.
We keep the pixel values in a big global array called |vv|.
The variables |m| and~|n| keep track of the current number of rows
and columns in use.

The |dd| table contains density values assumed for the input,
indexed by single-character codes.

@<Glob...@>=
@!vv:array [0..max_m,0..max_n] of real; {pixel darknesses, from 0.0 to 1.0}
@!m:integer; {rows |0..m+1| of |vv| should contain relevant data}
@!n:integer; {columns |0..n+1| of |vv| should contain relevant data}
@!dd:array[text_char] of real;

@ All input codes give zero density, except |"0"| to~|"9"| and
|"A"| to |"V"|.

@<Set init...@>=
for i←first_text_char to last_text_char do dd[chr(i)]←0.0;
for i←"0" to "9" do dd[chr(i)]←1.0-(i-"0")/31.0;
for i←"A" to "V" do dd[chr(i)]←1.0-(i-"A"+10)/31.0;

@ The process of inputting pixel values is quite simple. We terminate the
program if anomalous values of |m| and~|n| occur.  The values are duplicated
at the top, left, right, and bottom in order to provide ``padding'' that will
be convenient in the pixel transformation process.
@↑padding@>

@<Input a picture, or terminate the program@>=
read(m);@+if (m≤0)∨(m≥max_m) then finish;
read_ln(n);@+if (n≤0)∨(n≥max_n) then finish;
for i:=1 to m do
	begin for j:=1 to n do
		begin read(c); vv[i,j]:=dd[c];
		end;
	read_ln;
	end;
for j:=1 to n do
	begin vv[0,j]:=vv[1,j]; vv[m+1,j]:=vv[m,j];
	end;
for i:=0 to m+1 do
	begin vv[i,0]:=vv[i,1]; vv[i,n+1]:=vv[i,n];
	end

@ The code just written makes use of three temporary registers that must
be declared:

@<Glob...@>=
@!i,@!j:integer; {current row and column}
@!c:char; {character read from input}
@!cc:ASCII_code; {ASCII equivalent of |c|}
@* Approximating the darkness levels.
We shall use the following model for estimating the effect of a
given bit pattern: If a pixel is black, the darkness is 1.0; if it
is white but at least one of its four neighbors is black, the darkness
is |zeta|; if it is white and has four white neighbors, the darkness
is zero.

@d white=0 {code for a white pixel with all white neighbors}
@d gray=1 {code for a white pixel with 1, 2, 3, or 4 black neighbors}
@d black=2 {code for a black pixel}

@<Constants...@>=
@!zeta=0.2; {assumed darkness of white pixel with a black neighbor}

@ Each input pixel has a real-valued density between 0.0 and 1.0.
An input pixel corresponds to an $8\times8$ array of output pixels whose
darknesses are determined by the rule just mentioned. We shall call the
input pixels Big Pixels.

Let $v_{i,j}$ be the given density value of Big Pixel |(i,j)|. Then
$v_{i,j}$ will also be the density value of little pixel $(8i-4,8j-4)$,
and the densities of other little pixels are determined by linear
interpolation. More precisely, the density of little pixel $(8i-4+a,8j-4+b)$
when |0≤a≤8| and |0≤b≤8| is defined to be
$$(a/8)\bigl[(b/8)[v_{i,j},v_{i,j+1}],(b/8)[v_{i+1,j},v_{i+1,j+1}]\bigr].$$
(The densities are essentially constant near the boundaries, because of
the ``padding'' described earlier.)

There isn't room to store all the little pixels in memory at once, but
it suffices to keep buffers for about a dozen rows near the current area
being computed.

@<Glob...@>=
@!eight_n:0..big_n; {|8n|}
@!ii:integer; {the buffer holds rows |8ii-7| through |8ii+4|}
@!pure_row:array[0..big_n] of real; {densities in row |8ii+4|}
@!prev_pure_row:array[0..big_n] of real; {densities in row |8ii-4|}
@!delta_row:array[0..big_n] of real; {|pure_row-prev_pure_row/8|}
@!buffer:array[-2..9,0..big_n] of real; {densities in twelve current rows}
@!darkness:array[-3..9,0..big_n] of white..black; {darknesses in buffer rows}

@ The |interpolate_in| procedure computes the little pixels in |pure_row|
by looking at the Big Pixels in a given row of~|vv|.

@<Basic procedures@>=
procedure interpolate_in(@!ii:integer);
var @!acc,@!del:real; {current value and increment}
@!j:0..big_n; {runs through little pixels}
@!jj:0..max_n; {runs through Big Pixels}
@!bb:0..8; {distance from ``pure'' value}
begin jj←1; del←0; bb←5; j←1; acc←vv[ii,0];
repeat pure_row[j]←acc; incr(bb); incr(j);
if bb<8 then acc←acc+del
else	begin bb←0; acc←vv[ii,jj]; incr(jj); del←(vv[ii,jj]-acc)/8.0;
	end;
until j>eight_n;
pure_row[0]←0.0; pure_row[eight_n+1]←0.0;
end;

@ Here is a procedure that ``rolls'' the buffer down eight lines:

@<Basic procedures@>=
procedure roll;
var @!j:0..big_n;
@!i:2..9;
begin for i←6 to 9 do for j←0 to eight_n+1 do
	begin buffer[i-8,j]←buffer[i,j]; darkness[i-8,j]←darkness[i,j];
	end;
for j←0 to eight_n+1 do darkness[-3,j]←darkness[5,j];
for j←0 to eight_n+1 do prev_pure_row[j]←pure_row[j];
incr(ii); if ii≤m then interpolate_in(ii+1);
for j←0 to eight_n+1 do
	begin buffer[9,j]←pure_row[j]; darkness[9,j]←white;
	delta_row[j]←(pure_row[j]-prev_pure_row[j])/8.0;
	end;
for i←8 downto 2 do for j←0 to eight_n+1 do
	begin buffer[i,j]←buffer[i+1,j]-delta_row[j]; darkness[i,j]←white;
	end;
end;

@ It's tedious but not difficult to get everything started.
We put zeros above the top lines in the picture.

@<Initialize the buffers@>=
eight_n←8*n; ii←0; interpolate_in(1);
for i←6 to 9 do for j←0 to eight_n+1 do
	begin buffer[i,j]←pure_row[j]; darkness[i,j]←white;
	end;
for i←-2 to 5 do for j←0 to eight_n+1 do
	begin buffer[i,j]←0.0; darkness[i,j]←white;
	end;
for j←0 to eight_n+1 do darkness[-3,j]←white

@ It's easy to output the current darkness values. Here we output
eight consecutive rows.

@<Output the pixel values for the top eight rows of the buffer@>=
for i←-2 to 5 do
	begin write('row(',8*ii-5+i:1,'); cols('); cols_out←0;
	for j←1 to eight_n do if darkness[i,j]=black then
		begin if cols_out<15 then
			begin if cols_out>0 then write(',');
			incr(cols_out);
			end
		else	begin write_ln(','); write('               '); cols_out←1;
			end;
		write(j:1);
		end;
	write_ln(');')
	end

@ @<Glob...@>=
@!cols_out:0..15; {the number of columns output so far on this line}
@* Error diffusion.
Little pixels are divided into 64 classes, numbered from 0 to~63. 
We convert the pixel values by assigning them for class~0 first, then class~1,
etc. The error incurred at each step is distributed to the neighbors whose
class numbers are higher. This is done by means of precomputed tables
|class_row|, |class_col|, |start|, |del_i|, |del_j|, and |alpha| whose
function is easy to deduce from the following code:

@<Choose pixel values and diffuse the errors in the buffer@>=
for k←0 to 63 do
	begin i←class_row[k]; j←class_col[k];
	while j≤eight_n do
		begin @<Decide the color of pixel |[i,j]| and the resulting |err|@>;
		for l←start[k] to start[k+1]-1 do
			begin u←i+del_i[l]; v←j+del_j[l];
			buffer[u,v]←buffer[u,v]+err*alpha[l];
			end;
		j←j+8;
		end;
	end

@ @<Glob...@>=
@!class_row:array[0..63] of -2..8;
	{buffer row containing pixels of a given class}
@!class_col:array[0..63] of 1..8;
	{first column containing pixels of a given class}
@!class_number:array[-2..9,0..9] of 0..63; {number of a given position}
@!err:real; {error introduced at the current position}
@!err_black:real; {error introduced at the current position if black chosen}
@!black_diff:real; {difference between |err| and |err_black| for gray pixel}
@!l:0..256; {index into diffusion tables}
@!start:array[0..64] of 0..256;
	{first entry of diffusion table for a given class}
@!del_i,@!del_j:array[0..256] of -1..1; {neighboring location for diffusion}
@!alpha:array[0..256] of real; {constant of proportionality for diffusion}

@ Here we choose white or black, whichever minimizes the magnitude of the error.
The |gray| values of this pixel and its neighbors make this calculation slightly
tricky, as we must subtract |zeta| when a gray pixel is created and add |zeta|
when it is destroyed.

@<Decide the color of pixel |[i,j]| and the resulting |err|@>=
if darkness[i,j]=gray then
	begin err←buffer[i,j]-zeta; err_black←err-black_diff;
	end
else	begin err←buffer[i,j]; err_black←err-1.0;
	end;
if darkness[i-1,j]=white then err_black←err_black-zeta;
if darkness[i,j-1]=white then err_black←err_black-zeta;
if darkness[i,j+1]=white then err_black←err_black-zeta;
if darkness[i+1,j]=white then err_black←err_black-zeta;
if err_black+err>0 then
	begin err←err_black; darkness[i,j]←black;
	if darkness[i-1,j]=white then darkness[i-1,j]←gray;
	if darkness[i,j-1]=white then darkness[i,j-1]←gray;
	if darkness[i,j+1]=white then darkness[i,j+1]←gray;
	if darkness[i+1,j]=white then darkness[i+1,j]←gray;
	end

@ @<Set init...@>=
black_diff←1.0-2.0*zeta;
@* Computing the diffusion tables.
The tables for dot diffusion could be specified by a large number
of boring assignment statements, but it is more fun to compute them by a method
that shows some of the mysterious underlying structure.

@<Initialize the diffusion tables@>=
@<Initialize the class number matrix@>;
@<Compile ``instructions'' for the diffusion operations@>

@ The order of classes
used here is the order in which pixels might be blackened in a font
for halftones based on dots in a 45$↑\circ$ grid.

@<Basic procedures@>=
procedure store(@!i,@!j:integer); {establish new |class_row|, |class_col|}
begin if i<1 then i←i+8@+else if i>8 then i←i-8;
if j<1 then j←j+8@+else if j>8 then j←j-8;
class_number[i,j]←k; class_row[k]←i; class_col[k]←j; incr(k);
end;
@#
procedure store_eight(@!i,@!j:integer); {rotate and shift for eight classes}
begin store(i,j); store(i-4,j+4); store(5-j,i); store(1-j,i-4);@/
store(4+j,1-i); store(j,5-i); store(5-i,5-j); store(1-i,1-j);
end;

@ @<Initialize the class number matrix@>=
k←0; store_eight(7,2); store_eight(8,3); store_eight(8,2); store_eight(8,1);@/
store_eight(1,4); store_eight(1,3); store_eight(1,2); store_eight(2,3);@/
for i←1 to 8 do
	begin class_number[i,0]←class_number[i,8];
	class_number[i,9]←class_number[i,1];
	end;
for j←0 to 9 do
	begin class_number[-2,j]←class_number[6,j];
	class_number[-1,j]←class_number[7,j];
	class_number[0,j]←class_number[8,j];
	class_number[9,j]←class_number[1,j];
	end

@ The tricky part of this process is the fact that some values near the
bottom of the buffer aren't ready for processing until errors have been
diffused from the next bufferload. In such cases we go up eight rows
to process a value that has been held over.

@<Glob...@>=
@!hold:array[0..9,0..9] of boolean; {is this value too close to the bottom
	of the buffer to allow immediate processing?}

@ The ``compilation'' in this step simulates going through the diffusion
process the slow way, and records the actions it does (so that they
can all be done a high speed later).

@<Compile...@>=
for j←0 to 9 do hold[9,j]←true;
for i←0 to 8 do for j←0 to 9 do hold[i,j]←false;
l←0; k←0;
repeat i←class_row[k]; j←class_col[k]; w←0; start[k]←l;
for u←i-1 to i+1 do for v←j-1 to j+1 do
 if class_number[u,v]>k then
	begin del_i[l]←u-i; del_j[l]←v-j; incr(l);
	if u=i then w←w+2 {neighbors in the same row get weight 2}
	else if v=j then w←w+2 {neighbors in the same column get weight 2}
	else w←w+1; {diagonal neighbors get weight 1}
	end
 else	if hold[u,v] then hold[i,j]←true;
if hold[i,j] then class_row[k]←i-8;
@<Compute the |alpha| values for class |k|, given the total weight |w|@>;
incr(k);
until k=64;
start[64]←l

@ @<Compute the |alpha| values for class |k|, given the total weight |w|@>=
for ll←start[k] to l-1 do
	begin if del_i[ll]=0 then alpha[ll]←2.0/w
	else if del_j[ll]=0 then alpha[ll]←2.0/w
	else alpha[ll]←1.0/w;
	end

@ @<Glob...@>=
@!ll:0..256; {loop index}
@!u,@!v:integer; {neighbors of |i| and |j|}
@!w:integer; {the weighted number of high-class neighbors}
@!k:0..64; {the current class being considered}
@* The main program.
Now we're ready to put all the pieces together.

@<The main program@>=
@<Initialize the diffusion tables@>;
write_ln('input btone'); write_ln;
@<Input a picture...@>;
write_ln('width:=',8*n:1,';');
@<Initialize the buffers@>;
repeat @<Choose pixel values and diffuse the errors in the buffer@>;
if ii>0 then @<Output the pix...@>;
roll;
until ii>m;
write_ln('end.');
cleanup_and_terminate:
@* Index.
Here are the quantities declared and/or used in the program.
(The uses of single-letter variables aren't indexed.)